Multiscale implementation of infinite-swap replica exchange molecular dynamics.

Replica exchange molecular dynamics (REMD) is a popular method to accelerate conformational sampling of complex molecular systems. The idea is to run several replicas of the system in parallel at different temperatures that are swapped periodically. These swaps are typically attempted every few MD steps and accepted or rejected according to a Metropolis-Hastings criterion. This guarantees that the joint distribution of the composite system of replicas is the normalized sum of the symmetrized product of the canonical distributions of these replicas at the different temperatures. Here we propose a different implementation of REMD in which (i) the swaps obey a continuous-time Markov jump process implemented via Gillespie's stochastic simulation algorithm (SSA), which also samples exactly the aforementioned joint distribution and has the advantage of being rejection free, and (ii) this REMD-SSA is combined with the heterogeneous multiscale method to accelerate the rate of the swaps and reach the so-called infinite-swap limit that is known to optimize sampling efficiency. The method is easy to implement and can be trivially parallelized. Here we illustrate its accuracy and efficiency on the examples of alanine dipeptide in vacuum and C-terminal ß-hairpin of protein G in explicit solvent. In this latter example, our results indicate that the landscape of the protein is a triple funnel with two folded structures and one misfolded structure that are stabilized by H-bonds.

Kinetics of O2 Entry and Exit in Monomeric Sarcosine Oxidase via Markovian Milestoning Molecular Dynamics.

The flavoenzyme monomeric sarcosine oxidase (MSOX) catalyzes a complex set of reactions currently lacking a consensus mechanism. A key question that arises in weighing competing mechanistic models of MSOX function is to what extent ingress of O2 from the solvent (and its egress after an unsuccessful oxidation attempt) limits the overall catalytic rate. To address this question, we have applied to the MSOX/O2 system the relatively new simulation method of Markovian milestoning molecular dynamics simulations, which, as we recently showed [ Yu et al. J. Am. Chem. Soc. 2015 , 137 , 3041 ], accurately predicted the entry and exit kinetics of CO in myoglobin. We show that the mechanism of O2 entry and exit, in terms of which possible solvent-to-active-site channels contribute to the flow of O2, is sensitive to the presence of the substrate-mimicking competitive inhibitor 2-furoate in the substrate site. The second-order O2 entry rate constants were computed to be 8.1 × 10(6) and 3.1 × 10(6) M(-1) s(-1) for bound and apo MSOX, respectively, both of which moderately exceed the experimentally determined second-order rate constant of (2.83 ± 0.07) × 10(5) M(-1) s(-1) for flavin oxidation by O2 in MSOX. This suggests that the rate of flavin oxidation by O2 is likely not strongly limited by diffusion from the solvent to the active site. The first-order exit rate constants were computed to be 10(7) s(-1) and 7.2 × 10(6) s(-1) for the apo and bound states, respectively. The predicted faster entry and slower exit of O2 for the bound state indicate a longer residence time within MSOX, increasing the likelihood of collisions with the flavin isoalloxazine ring, a step required for reduction of molecular O2 and subsequent reoxidation of the flavin. This is also indirectly supported by previous experimental evidence favoring the so-called modified ping-pong mechanism, the distinguishing feature of which is an intermediate complex involving O2, the flavin, and the oxidized substrate simultaneously in the cavity. These findings demonstrate the utility of the Markovian milestoning approach in contributing new understanding of complicated enyzmatic function.

Resorcinol Crystallization from the Melt: A New Ambient Phase and New "Riddles".

Structures of the α and ß phases of resorcinol, a major commodity chemical in the pharmaceutical, agrichemical, and polymer industries, were the first polymorphic pair of molecular crystals solved by X-ray analysis. It was recently stated that "no additional phases can be found under atmospheric conditions" (Druzbicki, K. et al. J. Phys. Chem. B 2015, 119, 1681). Herein is described the growth and structure of a new ambient pressure phase, ε, through a combination of optical and X-ray crystallography and by computational crystal structure prediction algorithms. α-Resorcinol has long been a model for mechanistic crystal growth studies from both solution and vapor because prisms extended along the polar axis grow much faster in one direction than in the opposite direction. Research has focused on identifying the absolute sense of the fast direction-the so-called "resorcinol riddle"-with the aim of identifying how solvent controls crystal growth. Here, the growth velocity dissymmetry in the melt is analyzed for the ß phase. The ε phase only grows from the melt, concomitant with the ß phase, as polycrystalline, radially growing spherulites. If the radii are polar, then the sense of the polar axis is an essential feature of the form. Here, this determination is made for spherulites of ß resorcinol (ε, point symmetry 222, does not have a polar axis) with additives that stereoselectively modify growth velocities. Both ß and ε have the additional feature that individual radial lamellae may adopt helicoidal morphologies. We correlate the appearance of twisting in ß and ε with the symmetry of twist-inducing additives.

Structure modeling, synthesis and X-ray diffraction determination of an extra-large calixarene-based coordination cage and its application in drug delivery.

An extra-large octahedral coordination cage (CIAC-114) was designed and modeled based on Co4-p-tert-butylsulfonylcalix[4]arene (Co4-(SC4A-SO2)) subunits and 4,4',4''-(benzene-1,3,5-triyl-tris(benzene-4,1-diyl))tribenzoate (BBB) ligands via the isomorphic replacement approach built from an analogous cage structure with a smaller size. X-ray crystallography revealed that the crystals obtained through solvothermal synthesis exhibited the anticipated structure. Each CIAC-114 cage is assembled by six tetranuclear Co4-(SC4A-SO2) units as vertices, which bear a four-fold rotational symmetry, and eight tripodal BBB ligands as linkers, which hold a D3h symmetry. Among its analogues CIAC-114 has the largest overall peripheral diameter of 5.4 nm and an internal cavity of 2.7 nm. After treatment by supercritical CO2, a single crystal sample of CIAC-114 transformed into amorphous material with the retention of the cage skeleton, which demonstrated good adsorption properties towards a small drug molecule, ibuprofen (Ibu), evidenced by IR spectroscopy, (1)H NMR spectroscopy, N2 sorption analysis, and drug release experiments. The Ibu release experiment in phosphate buffered saline solution (pH = 7.4) revealed that CIAC-114 exhibited a slow drug release behavior.

Locating landmarks on high-dimensional free energy surfaces.

Coarse graining of complex systems possessing many degrees of freedom can often be a useful approach for analyzing and understanding key features of these systems in terms of just a few variables. The relevant energy landscape in a coarse-grained description is the free energy surface as a function of the coarse-grained variables, which, despite the dimensional reduction, can still be an object of high dimension. Consequently, navigating and exploring this high-dimensional free energy surface is a nontrivial task. In this paper, we use techniques from multiscale modeling, stochastic optimization, and machine learning to devise a strategy for locating minima and saddle points (termed "landmarks") on a high-dimensional free energy surface "on the fly" and without requiring prior knowledge of or an explicit form for the surface. In addition, we propose a compact graph representation of the landmarks and connections between them, and we show that the graph nodes can be subsequently analyzed and clustered based on key attributes that elucidate important properties of the system. Finally, we show that knowledge of landmark locations allows for the efficient determination of their relative free energies via enhanced sampling techniques.

Full kinetics of CO entry, internal diffusion, and exit in myoglobin from transition-path theory simulations.

We use Markovian milestoning molecular dynamics (MD) simulations on a tessellation of the collective variable space for CO localization in myoglobin to estimate the kinetics of entry, exit, and internal site-hopping. The tessellation is determined by analysis of the free-energy surface in that space using transition-path theory (TPT), which provides criteria for defining optimal milestones, allowing short, independent, cell-constrained MD simulations to provide properly weighted kinetic data. We coarse grain the resulting kinetic model at two levels: first, using crystallographically relevant internal cavities and their predicted interconnections and solvent portals; and second, as a three-state side-path scheme inspired by similar models developed from geminate recombination experiments. We show semiquantitative agreement with experiment on entry and exit rates and in the identification of the so-called "histidine gate" at position 64 through which ≈90% of flux between solvent and the distal pocket passes. We also show with six-dimensional calculations that the minimum free-energy pathway of escape through the histidine gate is a "knock-on" mechanism in which motion of the ligand and the gate are sequential and interdependent. In total, these results suggest that such TPT simulations are indeed a promising approach to overcome the practical time-scale limitations of MD to allow reliable estimation of transition mechanisms and rates among metastable states.

Microscopic mechanisms of equilibrium melting of a solid.

The melting of a solid, like other first-order phase transitions, exhibits an intrinsic time-scale disparity: The time spent by the system in metastable states is orders of magnitude longer than the transition times between the states. Using rare-event sampling techniques, we find that melting of representative solids-here, copper and aluminum-occurs via multiple, competing pathways involving the formation and migration of point defects or dislocations. Each path is characterized by multiple barrier-crossing events arising from multiple metastable states within the solid basin. At temperatures approaching superheating, melting becomes a single barrier-crossing process, and at the limit of superheating, the melting mechanism is driven by a vibrational instability. Our findings reveal the importance of nonlocal behavior, suggesting a revision of the perspective of classical nucleation theory.

Order-parameter-aided temperature-accelerated sampling for the exploration of crystal polymorphism and solid-liquid phase transitions.

The problem of predicting polymorphism in atomic and molecular crystals constitutes a significant challenge both experimentally and theoretically. From the theoretical viewpoint, polymorphism prediction falls into the general class of problems characterized by an underlying rough energy landscape, and consequently, free energy based enhanced sampling approaches can be brought to bear on the problem. In this paper, we build on a scheme previously introduced by two of the authors in which the lengths and angles of the supercell are targeted for enhanced sampling via temperature accelerated adiabatic free energy dynamics [T. Q. Yu and M. E. Tuckerman, Phys. Rev. Lett. 107, 015701 (2011)]. Here, that framework is expanded to include general order parameters that distinguish different crystalline arrangements as target collective variables for enhanced sampling. The resulting free energy surface, being of quite high dimension, is nontrivial to reconstruct, and we discuss one particular strategy for performing the free energy analysis. The method is applied to the study of polymorphism in xenon crystals at high pressure and temperature using the Steinhardt order parameters without and with the supercell included in the set of collective variables. The expected fcc and bcc structures are obtained, and when the supercell parameters are included as collective variables, we also find several new structures, including fcc states with hcp stacking faults. We also apply the new method to the solid-liquid phase transition in copper at 1300 K using the same Steinhardt order parameters. Our method is able to melt and refreeze the system repeatedly, and the free energy profile can be obtained with high efficiency.

Sampling saddle points on a free energy surface.

Many problems in biology, chemistry, and materials science require knowledge of saddle points on free energy surfaces. These saddle points act as transition states and are the bottlenecks for transitions of the system between different metastable states. For simple systems in which the free energy depends on a few variables, the free energy surface can be precomputed, and saddle points can then be found using existing techniques. For complex systems, where the free energy depends on many degrees of freedom, this is not feasible. In this paper, we develop an algorithm for finding the saddle points on a high-dimensional free energy surface "on-the-fly" without requiring a priori knowledge the free energy function itself. This is done by using the general strategy of the heterogeneous multi-scale method by applying a macro-scale solver, here the gentlest ascent dynamics algorithm, with the needed force and Hessian values computed on-the-fly using a micro-scale model such as molecular dynamics. The algorithm is capable of dealing with problems involving many coarse-grained variables. The utility of the algorithm is illustrated by studying the saddle points associated with (a) the isomerization transition of the alanine dipeptide using two coarse-grained variables, specifically the Ramachandran dihedral angles, and (b) the beta-hairpin structure of the alanine decamer using 20 coarse-grained variables, specifically the full set of Ramachandran angle pairs associated with each residue. For the alanine decamer, we obtain a detailed network showing the connectivity of the minima obtained and the saddle-point structures that connect them, which provides a way to visualize the gross features of the high-dimensional surface.

Temperature-accelerated method for exploring polymorphism in molecular crystals based on free energy.

The ability of certain organic molecules to form multiple crystal structures, known as polymorphism, has important ramifications for pharmaceuticals and high energy materials. Here, we introduce an efficient molecular dynamics method for rapidly identifying and thermodynamically ranking polymorphs. The new method employs high temperature and adiabatic decoupling to the simulation cell parameters in order to sample the Gibbs free energy of the polymorphs. Polymorphism in solid benzene is revisited, and a resolution to a long-standing controversy concerning the benzene II structure is proposed.

How to promote sluggish electrocyclization of 1,3,5-hexatrienes by captodative substitution.

Hexatriene electrocyclization, if not disfavored by its harsh reaction conditions, can be highly useful for the synthesis of complex organic molecules. Herein we developed a two-layer ONIOM method which could predict the activation free energy of hexatriene electrocyclization with an accuracy of about 1.0 kcal/mol. Using this carefully benchmarked method, we calculated the activation free energies for a variety of substituted hexatrienes. It was found that extraordinarily rapid electrocyclization could occur for certain patterns of captodative substituted hexatrienes, including 2-acceptor-3-donor hexatrienes, 2-acceptor-5-donor hexatrienes, and 3-acceptor-5-donor hexatrienes. The activation free energies for these systems could be up to 10 kcal/mol lower than that of the unsubstituted hexatriene, and therefore, their electrocyclization could proceed smoothly even at room temperature. The mechanism for the captodative effect on hexatriene electrocyclization could be understood by calculating the affinity between the donor and acceptor group in the reactant state and transition state of the reaction. If the affinity was stronger in the transition state, captodative substitution would produce an extra acceleration effect. It was shown that our theoretical results were in excellent agreement with the experimental data from the recent synthetic studies of hexatriene electrocyclizations. Thus, the theoretical tools developed in the present study could be used to predict not only how to accelerate the hexatriene electrocyclization via substituent manipulation but also under what conditions each particular electrocyclization could be accomplished in the real experiment.

Quantum-chemical predictions of redox potentials of organic anions in dimethyl sulfoxide and reevaluation of bond dissociation enthalpies measured by the electrochemical methods.

A first-principle theoretical protocol was developed that could predict the absolute pK(a) values of over 250 structurally unrelated compounds in DMSO with a precision of 1.4 pK(a) units. On this basis we developed the first theoretical protocol that could predict the standard redox potentials of over 250 structurally unrelated organic anions in DMSO with a precision of 0.11 V. Using the two new protocols we systematically reevaluated the bond dissociation enthalpies (BDEs) measured previously by the electrochemical methods. It was confirmed that for most compounds the empirical equation (BDE = 1.37 pK(HA) + 23.1E(o) + constant) was valid. The constant in this equation was determined to be 74.0 kcal/mol, compared to 73.3 kcal/mol previously reported. Nevertheless, for a few compounds the empirical equation could not be used because the solvation energy changed dramatically during the bond cleavage, which resulted from the extraordinary change of dipole moment during the reaction. In addition, we found 40 compounds (mostly oximes and amides) for which the experimental values were questionable by over 5 kcal/mol. Further analyses revealed that all these questionable BDEs could be explained by one of the three following reasons: (1) the experimental pK(a) value is questionable; (2) the experimental redox potential is questionable; (3) the solvent effect cannot be neglected. Thus, by developing practical theoretical methods and utilizing them to solve realistic problems, we hope to demonstrate that ab initio theoretical methods can now be developed to make not only reliable, but also useful, predictions for solution-phase organic chemistry.